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Abstract 
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Using methods of kinetic theory and liquid state theory we propose a description of the non- 



equilibrium behavior of molecular fluids which takes into account their microscopic structure and 
thermodynamic properties. The present work represents an alternative to the recent dynamic 
density functional theory which can only deal with colloidal fluids and is not apt to describe the 
hydrodynamic behavior of a molecular fluid. The method is based on a suitable modification 
of the Boltzmann transport equation for the phase space distribution and provides a detailed 

O ! 

description of the local structure of the fluid and of the transport coefficients under inhomogeneous 
conditions. Finally, we propose a practical scheme to solve numerically and efficiently the resulting 
kinetic equation by employing a discretization procedure analogous to the one used in the Lattice 
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£T) ■ Boltzmann method. 
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PACS (47.11.-j,47.61.-k,61.20.-p) 

I. INTRODUCTION 

Understanding the dynamics of fluids under non-equilibrium conditions is of capital im- 
portance for fundamental statistical physics as well as applied disciplines such as engineering, 
fluid mechanics, rheology, and physiology. In recent years there has been an upsurge of in- 
terest towards the study of transport phenomena in strongly inhomogeneous systems mostly 
motivated by important physical and technological applications such as microfluidics 1 , col- 
loids, oil recovery, lab-on-a-chip devices etc. These examples require the knowledge of struc- 
tural and dynamical fluid properties in the presence of restricted geometries and/or struc- 
tured substrates and of external gradients or time dependent fields. Typical time and length 
scales involved can be significantly shorter than those usually assumed in standard ther- 
modynamic and hydrodynamic theories. In order to go beyond these descriptions, several 
theoretical approaches have been developed which differ not only by the nature of systems 
under scrutiny (colloidal systems have different dynamical behaviors in comparison with 
simple fluids), but also from subjective factors such as the individual scientific background 
and personal taste. These methods include the dynamic density functional approach, the 
kinetic approach, mesoscopic methods and methods based on effective free energies. In the 
present paper we shall be concerned in some detail only with the first two. 

In the last thirty years, massive efforts have been devoted to develop techniques to study 
the properties of non-uniform interacting many particle systems, among these Density Func- 
tional theory (DFT) being perhaps the most versatile^. In DFT the single-particle density 
profile p(r) is obtained via minimization of the Grand potential functional. Functional 
derivatives of the Grand potential determine the multi-particle correlation functions and all 
the structural and thermodynamic properties for a system with arbitrary inhomogeneity. 
Recently, dynamical generalizations of these equilibrium methods have been applied to non- 
equilibrium problems such as diffusion, Stokes drift, polymeric fluids confined to cavities, 
etc. The dynamical density functional ( DDF) is apt to describe the relaxation of Brownian 
particles in a medium and can be applied in situations where the local velocity distribution 
is very close to the Maxwellian and the density field varies slowly in tinted. The method is 
based on the assumption of instantaneous equilibrium, i.e. that the correlation functions at 
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a given instant are identical to those of the same equilibrium system having the same equi- 
librium density profile. In the case of overdamped dynamics, which characterizes colloidal 
particles immersed in a solvent, the evolution is mainly governed by structural rearrange- 
ments so that such an interplay is evident and is justified to describe the evolution of the 
system in terms only of its density and density correlations. It is therefore not surprising 
that in DDF the velocity field is slaved to the density field and does not play any autonomous 
role. 

On the other hand, in molecular fluids the momenta of the particles are not damped 
by the interaction with the solvent so that the total momentum is conserved and the mo- 
mentum current must be treated on equal footing as the local density-^. Hydrodynamics 
describes the non-equilibrium state of a system by means of a set of local variables which are 
density, momentum and temperature of the fluid^. However, hydrodynamics does not apply 
to phenomena which are not slowly varying in space and one has to consider a finer level 
of description, such as one based on a suitable generalization of the Boltzmann transport 
equation for /(r, v, t), the phase space density distribution of particles with position r and 
velocity v at time t. The modeling of the interactions in the transport equation depends on 
the nature of the fluid and on the degree of accuracy required. The widely used Boltzmann 
collision operator has been studied for several kinds of interaction potentials and gives pre- 
dictions for the transport coefficients, but does not provide an accurate representation of 
the thermodynamics and structural properties of fluids^^ 1 ^. For simple applications one 
can even approximate further the interaction term by a linear relaxation term, as proposed 
by Bhatnagar, Gross and Krook (BGK) 14 . A more refined approximation is needed if one is 
interested in a dynamical description taking into account both the equilibrium properties of 
the fluid, such as the equation of state and the equilibrium density profile in an external force 
field, and the transport properties under inhomogeneous conditions. The Revised Enskog 
theory (RET), originally introduced by Van Beijeren and Ernst for hard sphere systems^, 
has the ability to describe the local fluid structure within a kinetic formalism. It can serve 
as a reference theory to study fluids with different types of interactions. However, even 
in the RET case a numerical solution is often too demanding in terms of computer speed 
and memory so that one has to resort to some simplifications. One of these has been pro- 
posed more than a decade ago by Dufty, Santos and Brey^ and consists in separating the 
slowly evolving part of /(r, v, i), associated with the five hydrodynamic modes, from the 
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fast non-hydrodynamic modes. 

In the present work, we will show how to obtain a workable numerical method based 
on this equation which includes information about the microscopic nature of the fluid and 
contains hydrodynamics as a limiting case. This paper is organized as follows: In section HT1 
we introduce the transport equation and the collision operator, we define the hydrodynamic 
fields and derive the balance equations. In section ITO1 we briefly derive, for the sake of com- 
parison, the DDF equation, in section [IV] we derive the equation of evolution for molecular 
fluids which differs from the DDF equation because it considers the evolution of the density, 
velocity and temperature fields altogether. In Sec. [V] we propose to use the Lattice Boltz- 
mann method as a strategy to solve numerically the transport equation. Finally, in Sec. IVII 
we summarize and draw some conclusions. 



II. KINETIC THEORY 



We consider a simple fluid, whose elementary constituents, the molecules, mutually in- 
teract via a pairwise, spherically symmetric potential U(r — r'). The statistical description 
of such a system is based on the exact BBGKY hierarchy of dynamical equations for the 
reduced distribution functions, whose first level is the following kinetic equation 



d t f(r, v, t) + v • V/(r, v, t) + • ^-/(r, v, t) = Q(f\r, v, t) + B(f\r, v, t) (1) 

where /(r, v, t) is the one particle phase space density distribution at time t and at the 
point (r, v), F(r) is an external velocity-independent force field, Q(f\r,v,t) represents the 
effect on the single particle distribution function of the interactions among fluid particles and 
B(f\r, v, t) is a coupling to an external agent, usually termed heat-bath. The interaction 
term is given by the following exact expression 

fi(r,v,t) = i-V v y dr' J dv7 2 (r,v,r',v / ,t)V r f/(|r-r / |) (2) 

involving the two particles distribution function, f 2 (r,v,r' ,v' ,t), which in turn depends on 
the three particle correlation function. However, in the simplest closure schemes fi(/) (and 
B(f)) can be expressed in terms only of f(r,v,t) so that ([Tj) becomes self-consistent and 
one can devise practical schemes of solution. 
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In order to avoid the difficult calculation of the two-particle distribution function, one 
can introduce an approximate truncation of the BBGKY hierarchy. That is, we assume the 
following factorization of f 2 

/ 2 (r, v, r', v', t) = /(r, v, t)f(r', v', t)g 2 (r, r', t\n) (3) 

where g 2 is the equilibrium static pair correlation function, a functional of the local density. 

One can further coarse grain the description and use the kinetic equation (CD) to derive, 
by the usual Chapman- Enskog method, the macroscopic hydrodynamic equations, where a 
limited set, (d+2), fields, namely density, momentum current and energy density are assumed 
to represent the state of the fluid. To achieve this goal and obtain an autonomous set of 
equations by relating the currents of the hydrodynamics fields to the gradients of these fields, 
one needs to find first the form of the distribution /(r, v,t) when the system is perturbed 
from equilibrium. Having the perturbed form of f(r, v,t), one can compute the transport 
coefficients and thus relate currents and fields, without assuming the phenomenological 
constitutive relations. The transport coefficients are usually studied under conditions of 
almost uniform density, so that the spatial dependence of /(r, v, t) can be taken into account 
only at linear order. However, under many circumstances the presence of density gradients 
strongly influences the thermodynamic and dynamical properties of a fluid so that it seems 
important to study the interplay between structure and dynamics. 

A. Hydrodynamic fields and balance equations 

A direct solution of eq. ([T]) in terms of the unknown /(r, v, t) is clearly beyond reach. We 
thus follow a different strategy, originally proposed by Dufty, Santos and Brey in a study 
of the hard-sphere system. Their approximation is tantamount to treat separately and 
accurately the part of contributing to the evolution equations for the hydrodynamic 
fields from the part relative to the non-hydrodynamic moments of the distribution, which 
can be safely approximated by a simple relaxation time ansatz. 

A convenient way to analyze eq. (TjQ) is to consider the equation for the lowest velocity 
moments of the distribution function, which correspond to the standard five hydrodynamic 
fields describing the slowly evolving part of /(r, v, i). We first introduce the average local 



5 



density 

n(r,t) = I dv/(r,v,t), (4) 



the average local fluid velocity 

u(r,t)= i— v / dw/(r,v,t) (5) 
n(r,i) y 

and the average local temperature 

where fee is the Boltzmann constant. We shall employ throughout the paper the Einstein 
summation convention that repeated indices are implicitly summed over and the notation di 
to indicate the partial derivative w.r.t. the i — th cartesian component of the vector r and 
d t to indicate the partial derivative w.r.t. time. 

A set of balance equations are obtained for density and momentum by multiplying both 
sides of eq. (CQ) by 1, and v, respectively, and integrating w.r.t. velocity 

d t n{r, t) + di{n{r, t) Ul (r, £)) = (7) 

and 

mn(r, t) [d tUj (r, t) + ^(r, t)d iUj (r, t)] + d t P t f\r, t) - F i (r)n(r, t) - Cf\r, t) = bf\r, t) (8) 

An analogous balance equation for the temperature field can be derived by multiplying 
(CQ) by m(v — u) 2 /2 and integrating w.r.t. v 

^n(r, t) [d t T(r, t) + Ul (r, t)^T(r, *)] + (r, t)d iUj (r, t) + (r, t) - (r, t) = b {2) (r, t) 

(9) 

To establish the momentum and temperature equations (jSJ) and ([9]) we introduced some 
quantities, which in general cannot be expressed in terms of the hydrodynamic fields. These 
are the kinetic components of the pressure tensor, indicated with the superscript "K", 

P l f\r,t)=m J dwf{r,w,t){v-u) t {v-u) j (10) 

and the kinetic components of the heat flux vector 

gf ) (r,t) = ^ J dv/(r,v,i)(v-u) 2 0;-^ (11) 



In addition, we denned two terms stemming from molecular interactions (recall that j Qdv 
and f Bdv = 0, because the number of particles is conserved) 



Cf ) (r,t) = m J dv{v-u)iQ{f\T,-v,t). (12) 



and 



C^(r,t) = jj dv(v-u) 2 n(f\r,v,t) (13) 
and two terms stemming from the coupling with the heat bath 

6 i (1) (r,t) = m J dv(v-u)iB(f\T,v,t). (14) 

and 

6«(r,t) = ^ J dv(v-n) 2 B(f\v,v,t) (15) 

Following ref.— one can relate the term with the excess part of the pressure tensor , 
py (r, t) , stemming from the interactions, 

cU(v,t) = -d j p£\v,t), (16) 

where the components of the excess part (over the ideal gas) of the pressure tensor can 
be computed using a method due to Irving and Kirkwood which gives the following exact 
expression 

PW(T,t) = -~ fdX [ dr 12 [ dvdv' Y 1 ? A(r+(l~A)r 12 , v, r-Ar 12) V, t) (17) 

2 Jo J J |r 12 | a|ri 2 | 

with r 12 = (r — r'). With such an identification, eq. (jSj) assumes the form of the standard 
macroscopic equation for momentum balance. In general, no analogous relation exists be- 
tween C^(r, t) and the excess component of the heat flux, and one should instead take into 
account also the transfer of potential energy in order to derive an expression in terms of 
macroscopic fluxes. 



III. COLLOIDAL DYNAMICS AND DENSITY FUNCTIONAL STRATEGY 

Let us go back to the non-equilibrium case. Equation ([T|) contains as limiting cases the 
Hamiltonian dynamics and the fully underdamped dynamics when the heat bath exerts a 
large friction on the particles. In the standard description of the behavior of colloids it is 
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assumed that the particles are subject to two kinds of forces from the surrounding fluid 1) 
a deterministic friction, proportional to their velocity and described by Stokes law and 2) a 
stochastic force, characterized by a white noise spectrum, resulting from the interactions with 
the solvent molecules. In formulas, each particle is subjected to the following combination 



of solvent forces — rwyv + ^2 / ymkBT^, where £ is a Gaussian white noise and 7 a friction 
coefficient. A fluctuation-dissipation relation has been assumed between the noise amplitude 
and the friction coefficient so that the steady state velocity distribution is Maxwellian — . The 
effect of these forces on the particle phase space distribution /(r, v, t) can be represented by 
the following term in eq. (fl]) 21 

S(-*)(/|r, v, t) = 7[~^ + ^ • v]/(r, v, t) (18) 

Using the specific form of the heat-bath term, eq. ffT8l) . we shall now recover the dynamic 
density functional equation. In this case b^\r,t) = — 771771 (r, t)ui(r, t). The reduction of 
eq. ([T]) to a DDF dynamics can be done systematically by employing a multiple time scale 
approach, as used by Marconi and Tarazonar^ 1 ^ 12 ^. The method eliminates systematically 
all high moments of the velocity distribution function except the zeroth moment, that is 
the number density, by employing an asymptotic expansion in the inverse friction parameter 
7 _1 . The result is an equation for the density field containing corrections to the standard 
non-linear diffusion equation usually employed in DDF. However, for the present scope, we 
shall employ an heuristic and much simpler method to contract the phase space description 
of eq. ([T]) into the DDF description. 

We start by considering an approximate solution of the transport equation of the form 

f DDF (r, v,t) = n(r,t)[l + V " U ^ t) ]0 o (v) (19) 

with o (v) = [2^y 3/2 exp(-f|^ and v\ = k B T/m. Substituting eq. (EE]) into (□) 

0o(v){[<9 t n(r,t) + di(n(r,t)ui(r,t)] + 
v- F 

-^-[d t (n(r, t)ui(r, t)) + v^dMr, t) -n(r, t) + jn(r, t)uAr, t)] + 

m 

(^T--S ij )[d i (u j (T,t)n(T,t))-— Ui (r,t)n(r,t)} =fi(r,v,t) (20) 

m 

Multiplying by 1 and by eq. (!20|) and integrating w.r.t. velocity, we obtain the continuity 
equation (j7]) and the following equation for the momentum current 

m^(n(r,t)u i (r,t)) + 7 mn(r,t) Mj (r,t) + ^T9 i n(r,t)(r,t)-F J (r)n(r,t)-Cf ) (r,^ = 0, (21) 

8 



The DDF equation is recovered if one takes 7 large and neglects the time derivative of the 
momentum current, n(r, t)u{(r, t)), viz. one neglects the time derivative in the l.h.s of eq. 
(12Tj) . Using the continuity equation ([7]) to eliminate the current one obtains the following 
diffusion equation 

d t n(r,t) = ^-diikBTdin^t) - Fi(r)n(r,t) - C^(r,t)} (22) 
In order to identify the full DDF equation one must specify the form of the collision term 

We begin by a simple ansatz for the interaction term ([2]) by assuming the following 
adiabatic approximation 

n DDF (r,v,t) = / dr> [ dy , f(r,y,t)f(r',V,t)g 2 (r,r f ,t\n)^-U(\r-r'\) (23) 

m ovi J J OTi 

where we have approximated the two-particle correlations by those of an equilibrium system 
having the same density profile as the system at time t. Next, we define the molecular field 

as 

Ft° l \r,t) = - J dv'n{v\t)g 2 {v,v\t\n)^-U{\r-v'\) (24) 
so that the interaction term becomes 

n DDF {r,v,t) = --^-/(r,v,t)F/ mo °(r,t) (25) 
By computing the corresponding C 4 from (112j) the result is 



C\ l \v 1 t)=n{v,t)Ft° l \v,t). (26) 



It is straightforward to show that expression ff26l) can be recast as 

where AJ 7 is the free energy excess over the ideal ga-s^ . 

Notice that, in the case treated above, both the noise and the friction are externally 
imposed by the presence of the solvent and such an heat bath term breaks the translational 
invariance of the system, since the heat bath is assumed to be at rest. 

Unfortunately, as shown in ref.— , the method cannot be extended to systems with small or 
vanishing friction proportional to the velocity and therefore cannot be extended to describe 



dense liquids, which are instead characterized by an internal friction proportional to velocity 
gradients and not to the velocity itself. The slaving of the momentum current and of the 
energy current to the density, which allows the reduction of the complex dynamics eq. (pQ) 
to eq. (|22p is at work only in the high friction regime. Hence, a different approach is needed 
to treat molecular fluids. 

IV. MOLECULAR FLUIDS AND MODIFIED BOLTZMANN EQUATION 
STRATEGY 

The dynamics of molecular fluids is described by Newton's equation of motion, as opposed 
to colloidal fluids where the presence of the solvent corresponds to overdamped Brownian 
dynamics encapsulated in the DDF equation. Two important properties of molecular fluids 
must be preserved in any theoretical representation. These are the Galilei invariance, with 
respect to the reference frames moving relative to each other at constant velocity, and 
the momentum conservation. These are relevant features of hydrodynamics and are not 
displayed by the DDF method, where the absolute velocities are damped, thus privileging 
the "solvent reference frame" . In the case of molecular fluids, the total momentum must be 
conserved, so that a heat bath such as that used in the previous section is not acceptable. 
In a molecular fluid the dissipation is determined by other mechanisms such as, among the 
others, the internal friction proportional to velocity gradients (viscosity), and not to velocity 
itself, and to temperature gradients (thermal conductivity). We shall use this information in 
order to make approximations not violating these conservation laws, which are at the origin 
of the low frequency hydrodynamic modes. 

We start by splitting Q into two contributions, the first taking into account accurately 
the effect of the interaction term on the hydrodynamic variables, so that a correct thermo- 
dynamic and structure of the fluid is achieved, and the second describing in lesser detail the 
evolution of the non-hydrodynamic modes. 

Based on the form of eqs.(0|, ([5]) and ([6]), we write 



n(/|r,v,t) = ^((v-^C^^j^^-llC^Mjj+ffii/lr.v.t)) (28) 



2 



with 




(29) 
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Equation (|28|) represents a separation of the interaction term fl(f) into a direct contribu- 
tion to the hydro dynamic modes, plus a correction, 5Q, which is assumed not to act directly 
on these modes. It is easy to see that the following approximation for the term 5Q fulfills 
such a condition 

*n(/|r, v, t) ~ VQ [ kBT ^ t) ^ + ± ■ (v - u(r, t))]f(r, v, t) (30) 

In fact, by direct inspection one sees that the form of fi(f) eq. fTSg]) together with eq. ffHU]) 
reproduces eqs. ( fT2l) and (TIB"]) , and that <50 does not contribute to these two integrals. The 
term <5f2 is assumed to be linear with respect to / and is approximated by an heat-bath 
operator. At variance with the heat bath operator introduced in the colloidal problem, its 
effect is to induce a matching of the velocity of the particles to the local value of the average 
velocity of the fluid and a matching of the temperature of the particles to the local value 
of the average temperature. The coefficient z/ is a phenomenological adjustable parameter 
chosen to reproduce the value of the viscosity of the fluid. The simplified and approximated 
version of ([I]) represented by eq. (f28l) with eq. ff30|) reproduces by construction equations 
([7]), ([8]) and (Q. We comment that the form (130]) preserves the Galilei invariance of the 
fluid. This heat bath "co-moving" with the fluid is of course very different from the one 
introduced in the study of colloids. Physically it can be thought to represent the effect of 
the fast modes on the slow modes, and can be seen as an intrinsic heat bath, as opposed to 
the extrinsic heat bath employed in the study of colloids. 

We remark that the invariance of the heat bath with respect to the choice of the reference 
frame is the same requirement which has lead to the formulation of the so-called Dissipative 
Particle Dynamics (DPD) 2 - 1 ^. There on a phenomenological basis one introduces pairwise 
dissipative and random forces, a "DPD thermostat" which locally conserves the momen- 
tum and leads the emergence of hydrodynamic flow effects on the macroscopic scaled. The 
present method instead of introducing a pairwise friction proportional to the velocity dif- 
ference of two colliding particles, employs a frictional force proportional to the difference 
between the individual particle velocity and the average fluid velocity u. 

It is useful to represent f(r, v, t) as the sum of a local thermodynamic equilibrium state, 
//oc(r, v, i) = ra(r, t)0M(r, v, t), plus a contribution representing the deviation from such a 
state 

/(r,v,t) = / loc (r,v,t) + 5/(r,v,t) (31) 
11 



The function /z oc (r, v, t) alone fully determines the values of the hydrodynamic moments fl4]), 
([5]) and §6§, whereas Sf(r, v, t) does not contribute. On the other hand, 5f does contribute 
to the heat flux and to the viscous pressure tensor as shown below. The interplay between 8f 
and the hydrodynamic fields occurs both via the terms and C^ 2 \ which are functionals 
of Sf and fi oc and via kinetic components of the pressure tensor and heat flux vector. 

It is relevant to remark that 5Q vanishes when / = fi oc . The mathematical properties 
of the differential operator featuring in eq. fl30l) are well-known, in particular displaying a 
non-positive spectrum with discrete eigenvalues 0, —Uq, — 2i/ , .... The local Maxwellian (j29|) 
represents the eigenfunction associated with the null eigenvalue, whereas the higher-order 
eigenf unctions are associated with the non hydrodynamic modes. However, for practical pur- 
poses it is more convenient to simplify further this term and resort to a drastic assumption, 
by replacing the heat bath term (1301 by a BGK-like relaxation term 

6Q(f\r,v,t) = -i/o(/(r,v,t)-n(r,t)^(r,v,t))= -v 6f(r, v, t) (32) 

which keeps the relevant properties of (13"U1) . at the price of assuming a single relaxation time, 
Vq 1 , for all non-hydro dynamic modes. 

Here, vq is a phenomenological parameter, representing a collision frequency, chosen as to 
reproduce the viscosity of the fluid. The above approximation replaces by a much simpler 
relaxation the complicated effects of the interactions among the non-hydro dynamic moments. 
The action of 5Q is to induce a rather fast relaxation of the solution towards a state of local 
equilibrium^. Once such a state has been reached the system evolves towards the steady 
state via equilibration of different regions through exchange of hydrodynamic fluxes. 

In addition, since also the term enforcing the local equilibrium condition (130]) contains as 
local parameters u(r, t) and T(r,t), the system has the correct long wavelength properties 
required by hydrodynamics. For future reference we rewrite the combination of eq.([T]), (128]) 
and 

d t f(v, v, t) + v ■ V/(r, v, t) + ^ • — f(r, v, t) - 

m av 

The advantage of using eq. (1331 instead of the apparently equivalent set of coupled hy- 
drodynamic equations is the following. Even in the non-interacting case the moments of 
the distribution function are coupled and to interrupt the hierarchy one needs a truncation 
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scheme. This truncation can be avoided by working directly with eq. (1331) . where one can 
apply the powerful technique of the Lattice Boltzmann equation (LBE). In the LBE one 
discretizes the velocity and the coordinates of the particles using a finite grid and computes 
directly the distribution function /(r, v, £). This program would had been numerically too 
demanding if we had to use eq. ([T]) because of the large number of integrals involved in the 
evaluation of the interaction term. 

Eq. fl33l) treats the hydrodynamic moments in a privileged fashion and describes a rapid 
equilibration of the system toward the local equilibrium state. The remaining stage is 
described by hydrodynamics, that is by mass, momentum and energy transport on larger 
scales. 



A. Approximate solution, kinetic contribution to the viscosity and heat conduc- 
tion 

Although this procedure is not necessary in numerical work, in order to gain some insight, 
in the present section we shall determine 5f perturbatively starting from fi oc , the local 
equilibrium state. We substitute / ~ j\ oc in the l.h.s. of (133!) and obtain an equation for Sf 
in terms of derivatives of the hydrodynamic fields. 

The details of the calculation can be found in standard textbooks 2 - and hence will be 
skipped here. The substitution of fi oc in the left hand side of eq. ( 1331) gives 

Fk ^(1)1 m(v k - u fe ) 



k B T 



| [d t n + di(nui)] + [nd t u k + nuidi(nu k ) + d k {nT) - — - C^] ■ 

n r _ m _ m 2 mr , 2 „ f , n ,m(v — u) 2 , n „ , (v — u) 2 5 X/ . 

+ — \d t T + uAT + -TdiUi C (2) l (—, - 3) + -diTim^— '- ) (v % - m) 

2T l 3 3n JV k B T ' T v 2k B T 2 n ' 

nr. w , (v — u) 2 \ 1 , / \ ^/( r ) v ,^) / nl , 
+m— — [{Vi -Ui){v k - u k ) dijjdiUk \(p M (r,v,t) = (34) 

K B 1 \ O J J h>Q 

Since Sf does not contain terms proportional to the first three terms of the above equation, 
because it must have vanishing hydrodynamic moments, we must impose that the first three 
terms in the l.h.s. vanish. These are the so-called solvability conditions and are precisely 
the balance equations ©, (JHJ) and (jUj) at the Euler level, i.e. without kinetic viscosity and 
heat conduction contributions. 
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We thus obtain the following explicit representation of 5f 



if = -^ J F i ^( r .-.')[(T^y-i)(«<-«-)^.') 

/ (v — u) 2 \ ] 

+ m(Ju i -«i)(t> fc -'Uj fc ) 5ijjdiU k {v,t) . (35) 



i/ T 

^ 2 

3 

The actual solution 5f, which can be determined by solving numerically eq. fl33l) contains 
higher-order terms, but we shall not try to go beyond the approximation (|35|) in the present 
paper. Instead, we determine the kinetic contribution to the transport coefficients. We first 
compute the heat flux vector by substituting 5f in eq. (|TTT) 



g 4 (r, t) = ~— n (r, t)k 2 B T(v, t)d<T(v, t) (36) 
and secondly we compute the components of the pressure tensor ffTUl) 

F#°(r,f) = fc fl T(r,t)n(r,t)(yy - t} ((<9^(r, t) + d jUi (r, t)) - 2 -d k u k {r,t)5^ 

(37) 

By comparing ( 1361) with the macroscopic expression g« = —XdiT we obtain the kinetic 
contribution to the heat conductivity 

A<*> = |— nfc|T (38) 

By comparing ( |37|) with the macroscopic definition of pressure tensor, 

P*>«\r, t) = k B Tn6 tJ - ^(d^ + d jUi ) + (r/f > - ^ K) )d k u k 5^ (39) 

where rjb is the bulk viscosity coefficient, we find the kinetic contribution to the shear viscosity 
coefficient 

= ^ (40) 

and 

vl K) = (41) 

It is important to stress that in the Euler approximation it is not possible to have a stationary 
solution because there is no dissipation mechanism to release the energy injected by an 
external force, so that viscosity and heat conduction are necessary. Finally, in order to 
make contact with the literature one can fix the free parameter of the theory, u , and 



choose Uq = na 2 '/m, with a the diameter of the equivalent hard sphere, and find 



(K) _ k B Tm 2 
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B. Short range repulsive potentials 



To proceed further, one must solve eq. (1331) and obtain explicit expressions for the 
thermodynamic fields and the transport coefficients. Thus, we need a specific form of the 
interaction potential U(r,r') and consequently of Q(r, v, t). We shall relate the microscopic 
details to the transport coefficients so that we need to compute the quantities Cj; and C^. 

The prototypical short range repulsive potential is represented by the hard sphere po- 
tential for which one has to consider a special treatment of the interaction term in order to 
obtain an accurate representation of the excess quantities. In particular, such an interaction 
can be treated as a collision process and the collisions as an uncorrelated binary sequence. 
A first approximation of fl[/](ri, v, t) is given by the "Stosszhal ansatz" which renders ([T]) 
an equation involving only the single-particle distribution decoupled from higher-order dis- 
tribution functions 

Q B [f](T,v,t) = a 2 J dv 2 J cffi0(6-v 12 )(6- v 12 ) x 

[/(r, v', t)f(v, v' 2 , t) - f(v, v, t)f(r, v 2 , t)} (42) 

where v' and v 2 are scattered velocities given by v' = v — (6 • v 12 )6 and v' 2 = v 2 + (6 • v 12 )6 
with V12 = v — v 2 . 

At higher densities, however, the "Stosszhal ansatz" fails to describe both the structural 
and relaxational properties of the fluid because collision sequences become highly correlated. 
To include these sequences and extend the transport equation to higher densities in the 
Seventies van Beijeren and Ernst have developed the revised Enskog theory (RET), taking 
into account the effects of ternary and higher-order collisions and the difference in positions of 
two hard-spheres at collision. Such a feature allows the instantaneous transfer of momentum 
and energy during a collision. In particular this transport mechanism gives rise to non-ideal 
gas contributions to the pressure and to the heat flux, which were neglected in Boltzmann's 
treatment of collisions. The RET collision operator takes the form 

CW/](r,v,t) = a 2 J dv 2 J cffi6(6 ■ v 12 )(6 ■ v 12 ) x 

\g 2 {v, r - 66, t\n)f(r, V, t)/(r - fifi, v' 2 , t) 

- g 2 (r, r + fifi, t\n)f(r, v, t)f(r + 66, v 2 , t)} (43) 

In order to obtain an explicit representation of the excess pressure tensor and heat flux we 
must compute, as prescribed by f|T2|) and ffTSl . the integrals of ^^^(r, v, t) times Gf\v) = 
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m(v — u(r, t))i and G (2) (v) = ~(v - u(r, t)) 2 . After replacing (v, v 2 ,fi) — ¥ ( v '> v 2; — ")> the 
RET explicit form of eqs. (I12II13P reads 

C®{r,t) = ^Jdvjdv 2 J dfi9(fi- v 12 )(fi- v 12 )[Gf (V) -Gf (v)] x 
<? 2 (r, r + fifi, t|n)/(r, v, t)f(r + fifi, v 2 , t) - ^(r, r - fifi, *|n)/(r, v 2 , f)/(r - fifi, v, t)ju) 

In the case of the hard-sphere fluid, where there is no contribution to the internal energy 
stemming from the pair potential, one finds the simple relation 



(45) 



where the first term represents the divergence of the heat flux and the second the contribution 
due to the viscous heating. 

At this stage we simplify drastically the calculation of the interaction term and neglect the 
contribution 6f(r, v, t) to the integrals (jSJ). We assume /(r, v, t) = n(r, i)0 M (r, v, t), which 
depends on the three hydrodynamic fields n, u, T. Substituting now this approximation into 
( |44l) and expanding to first order <pM (r ± fi, v, t) about u(r, t) and T(r, t) 



>A1 



(r ±fifi,v,t) 



m 



2nk B T { r,t) 



3/2 



exp 



m(v — u(r, t))' 
2k B T{r,t) 



(1 + 



m[v - u(r, t)] ■ [u(r ± fifi) - u(r)] 
k B T(r,t) 

f m[v-u(r,t)]-[v-u(r, t))-3^T(r,t)] [T(r ± ^ f) _ ^ t)] + 



2fc B T 2 (r,t) 

we arrive, after some lengthy algebra, to the result 



C^\r,t) = ~k B T{r,t)n{r,t)a 2 I dfifi^ 2 (r, r + fifi, t\n)n(r + fifi, t) 



T{r + W,t)-T(r,t) 



nk B T(r, t)/m 



fi • [u(r + fifi, t) - u(r, t)] + 



2T(r,t) 



(46) 



and 



C (2) (r, t) = k B T(r, t)n{r, t)a 2 J c% 2 (r, r + fifi, t\n)n{r + fifi, t) 
fi - [u(r + fifi,t) -u(r,t)] 



1 / k B T(r, t) T(r + fifi, t) - T(r, t) ' 
7r V m T(r,t) 



(47) 
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A comment is in order. In the hard-sphere fluid, momentum and energy fluxes can be 
transferred instantaneously even when the velocity distribution function has a Maxwellian 
form, provided it peaks at the local hydrodynamic velocity with a spread determined by the 
local temperature. As a consequence, we obtain a contribution not only to the pressure but 
also to the transport coefficients even within a Maxwellian approximation to the distribution 
function. This result is at variance with the corresponding result in the case of the Boltzmann 
equation where the collision term does not contribute neither to the pressure nor to the 
transport coefficients when the distribution is Maxwellian. 



C. Collisional contribution to the transport coefficients in the Hard-Sphere fluid 

We now apply our method to compute the fluid transport coefficients in bulk conditions. 

Shear viscosity. We assume both the density and the temperature to be uniform through- 
out the system, i.e. n(r, t) = hq and T(r, t) = Tq whereas the velocity profiles varies along 
a direction normal to the stream lines 

u = (0, u y {x, 0),0) (48) 

We then use the equations 
and 

r)p( C ) Pi 1 it 

c? = - 9J t = ^ («» 

in order to establish a relation between the macroscopic 77 coefficient and the microscopic 
level represented by . 

Using eq. ( H6i) and expanding to second order in a we obtain 

Cf (r, t) = mv *a%(cr)nZ-±-^ J d™»a& = mv T a% (<7)n »-L^il (51) 

where (72(c) is the radial distribution function at contact. By comparing ( 15T1) and ( l50i) we 
obtain 

rf c ^ = — mt>rV / 7r°" 4 fl , 2( cr )^o = 0.266a/ k B TTrmg 2 (o-)nla 4 (52) 

which is the viscosity found by Longuet-Higgins and Pople^ the case of hard-spheres, 

(c) 



whereas Enskog formula^ at high density gives ry E ^ skog = 0.337v / ^B7Vm(yf 2 (cr)nQ(T 4 . 
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Bulk viscosity. Similarly if one assumes u = (u x (x, 0),0,0), recalls cP(r,t) = and 
uses the macroscopic relation 

dP xx . 4 d 2 u x 
can extract the bulk viscosity. Expanding C x to second order in a we find 

4.\ 2 4 I \ 2 1 <9 2 ^a; f ,4^4 4 / \ 2 1 <9 2 ^x 47T 

Q j M) = "™t^ 92 ^' n ° yfUvT dx 2 J x = mVTa g2 ^' n °^/n~dx ir ~5 ^ ' 

That is r^ c) = §r/ c ). 

Thermal conductivity. The thermal conductivity can be computed assuming that both 
the density and the velocity are uniform throughout the system, i.e. n(r, t) = n and 
u(r, t) = 0. In this case the heat flux is proportional to the temperature gradient: q 1 - ^ (r, t) = 
-A^VT(r,t) and C^(r,t) = -V • q( c )(r,t). We assume a temperature profile varying 
only along the x-direction T(x, t) and use eq. (l4"Tl) . so that after some simple calculations and 
expanding to second order in a we obtain 

C^(r,t) = mv 3 T g 2 (a)ny±^~V 2 T(r,t) [ d d ~'aa 2 x = m, T AWn^^V 2 T(r, t) 

•y 7T TflVrp Z J Tfl o 

(54) 

The heat conductivity HS excess is 

= - \/ kBTTimg2{o-)nia 4: — = O.Q66^kBTTrmg 2 (a)n 2 ] o- A — (55) 
3 m m 

while Enskog formula at high density gives X^nskog = ^■^^\ / kBT7rmg 2 (a)nla 4 ^- 



D. Beyond hard core dynamics 

In the case of more realistic potentials, such as the Lennard- Jones interaction, a satisfac- 
tory theory for dense systems is still missing. The situation has been reviewed recently by 
Stell and coworkers^. They compared different theories of transport in the Lennard- Jones 
fluid. These include a kinetic variational theory (KVT)^, and a stochastic approach put 
forward in its original form by Leegwater— and later reformulated by Polewczak and Stell^ 2 . 
which introduces a random distribution of diameters of hard-spheres in the collision term. 

The KVT applies to Lennard- Jones (LJ) fluids and is obtained by adding a hard-sphere 
core to the attractive tail of the LJ potential. The transport coefficients exhibit Enskog- 
like forms, but the radial distribution function bears explicit dependence on the LJ tail as 
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well as on the hard-sphere core. The hard sphere diameter is determined according to the 
well-known WCA method used in equilibrium statistical mechanics to mimic the LJ fluid 33 . 

We shall not try to apply these theories in the present paper, a task beyond our scope, 
although the stochastic method could perhaps be implemented in our scheme. We shall only 
comment that the difficulty for continuous interactions is that there are not instantaneous 
collisions followed by free streaming trajectories as in the case of hard spheres. Instead, the 
collisions have a finite duration and a second collision event can take place before the first 
one is completed. A possible workaround could be to consider the motion of each particle 
as the combination of continuous momentum changes and almost instantaneous collisions. 
These small momentum changes represent the effect of weak attractive forces exerted by the 
surrounding molecules and can be assimilated to the random white noise forces occurring 
in the Brownian motion. In the literature, Rice and Allnatt combined the hard-core and 
the soft fluctuating potential model^ 1 ^. They represented the collisional term as a sum of 
the Enskog-Boltzmann collision integral, describing the large momentum binary exchanges 
plus a Fokker-Planck collision operator, supplemented by an average force, to model the 
small momentum exchanges due to the attractive forces, and assumed that these processes 
proceed without mutually interfering. 

This separation should be feasible, since in many simple fluids the pair potential can be 
decomposed into the sum of a short range harsh repulsive potential contribution, plus a an 
attractive long range tail. The fine details of such a decomposition depend on the accuracy of 
the approximation, the most popular methods being the Weeks- Chandler- Andersen (WCA) 
and the Barker-Henderson 3 ^ ones. It is therefore possible to treat the short range part with 
the technique employed in the case of the hard spheres and the tail of the interaction in 
the adiabatic approximation illustrated in section IHIl On the other hand, there are other 
possibilities. One of them consists in approximating the full continuous potential by a 
stepwise potential for which one can derive a proper extension of the RET equation, but the 
method seems exceedingly cumbersome to be of practical relevance. 

In conclusion, as far as we are interested in studying systems of particles with realistic 
potentials we can treat the attractive tails in a mean field (Vlasov) fashion and the repulsion 
as if determined by a hard-spheres with a diameter suitably chosen. Since and C^> 
do not have dependence on velocity and temperature gradients, they cannot contribute to 
the transport coefficients. To be more explicit, by projecting the interaction Q over the 
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hydrodynamic space, the resulting "matrix elements" do not depend on the fields u and 
T, i.e. the Vlasov approximation only accounts for the distortion of the density n from its 
equilibrium value. Thus, the transport coefficients are only related to the BGK term and 
are purely kinetic within this approximation. 

V. SOLUTION BY LATTICE BOLTZMANN STRATEGY 

The Lattice Boltzmann represents an alternative to methods which involve the direct 
solutions of coupled hydrodynamic equations of the type (j7|), (JS]) and fl5J). These equations 
represent projections of the evolution, eq. fl33l . and are not closed, since even the kinetic 
contributions of the pressure tensor and heat flux vector are not expressible solely in terms 
of the n, u and T. The necessary information to determine these contributions is contained 
in the higher moments of the distribution function. One can avoid the difficulty of computing 
these moments by using a macroscopic point of view which assumes constitutive relations, 
i.e. introduces phenomenological linear relations between and q and the velocity and 
temperature, respectively. Alternatively, one is faced with the problem of closing a dynamical 
infinite hierarchy of equations. We instead propose to solve eq. (1331) directly, with the only 
approximation stemming from the discretization procedure. 

The LB method first evolved as empirical extension of lattice gas automata and found 
wide application as a standard simulation tool in computational fluid-dynamics? 7 . In sub- 
sequent years, the LB method was found to be a systematic procedure to solve numerically 
a kinetic equation in velocity space^^ 1 ^. Having in mind the simulation of condensed sys- 
tems in micro/mesoscopic conditions, the critical parameter governing the departure from 
equilibrium is the Knudsen number, e, being the ratio between the mean free path and the 
representative length scale. On the nanoscale, e can be arbitrarily large and the numerical 
method should offer great flexibility with respect to external conditions and the Knudsen 
value. 

As previously shown, a direct solution of eq. ([1]) is numerically very demanding in the 
case of a realistic form of the collision term Q(f) due to the large number of integrations over 
the phase-space variables. In this section we focus on the general strategy to solve the kinetic 
equation (I3"3l in the adiabatic approximation, i.e. for = and the system is rigorously 
isothermal. However, extending these concepts to the presence of heat transport is rather 
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straightforward, along the lines discussed below. Moreover, we assume that the functional 
g2(r,r',t\n) is completely known by some level of theory or from atomistic simulations, 
such as Molecular Dynamics. In the study of fluids under non-equilibrium conditions, the 
proposed scheme has a strategic advantage over conventional atomistic simulations schemes 
since it exploits the pre-averaged nature of the kinetic description, avoiding the need to 
average observables over different realizations of the noise. 

The starting point to derive a consistent numerical scheme is the representation for both 
the distribution function and the collisional terms over a finite set of Hermite polynomials. 
The distribution f(r,v,t) is approximated by 




with u(v) = (2Trv^)~ 3 ^ 2 e~ v2 ^ 2v t . To keep the notation compact, in the following we use the 
convention that the product of tensors implicitly indicates the sum over all permutations of 
the tensorial indices. By construction, the complete and truncated distributions have the 
same coefficients up to Hermite order K, i.e. 

J c/v/i (/) (v)/(r,v,t) = J c/v/i w (v)/(r, v,t) for I < K (57) 

Being the distribution coefficients 0^ a combination of the distribution moments up to I 
order, the full and truncated distributions share the same moments up to K order. 

Similarly, the collision operator is replaced by a truncated representation over the 
same Hermite set. We distinguish between the non-hydro dynamic relaxation term, taken 
here as the BGK relaxation fl32|) . and the collisional term, renamed as /C(r,v, t) = 
±CP(r,v,t)d Vi $ M (r,v,t) = _^cf 1) (r ) v,t)«9, ; J e «(r,v,t). Alternatively, the straight 
heat-bath operator ( |30l) or momentum-preserving version ( l30l) can be chosen and treated 
along similar lines^ 1 ^. In this way, the colloidal/density functional dynamics or a different 
microscopic dynamics for the non-hydrodynamic modes can be selected. The corresponding 
truncated representations are 




and 

7C(r, v, t) = u(v) lLx (0 (r, t)h® (v) (59) 

1=0 V T AL - 
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Again, each collisional term shares the same K moments with the original collisional coun- 
terpart. 

The next step is to employ Gauss-Hermite quadratures to evaluate hydrodynamic and 
collisional moments. Recognizing that /(r, v, t)h^(v)/u(v) = p(r, v, t) is a polynomial in v 
of degree < 2K, the moments can be evaluated exactly with quadratures of order 2G > K, 
since 

w (r,t) = yrfv/(r,v,t)/i (z) (v) = J dvu(v)p(r,v,t) 

G G 
p=0 p=0 

where f p (r, t) = w p f(r, c p , t)/u(c p ), the vectors c p are a set of quadratures nodes and w p the 
associated quadrature weights. Similarly, the collisional moments are computed as 

G 

p=0 
G 

X«(r,t) = ]T/C p (r,t)/^(c p ) (62) 

p=0 

In summary, with the truncated Hermite representation and Gauss-Hermite quadra- 
tures, the distribution is replaced by an array of Q populations, /(r, v,t) — > f p (r,t) = 
w p f(r,c p ,t)/ui(c p ) and similarly for the collisional term, <5f2(r,v,t) — > 5Q p (r,t) = 
w p 5Q(r,c p ,t)/u}(c p ) and /C(r, v,t) ^ /C p (r, t) = w p JC(r, c p , t)/u(c p ). The nodes c p are chosen 
as vectors connecting neighboring mesh points r on a lattice, mirroring the hop of parti- 
cles between mesh points, generally augmented by a null vector c accounting for particles 
at rest. The specific form of the lattice velocities and weights depends on the order of 
accuracy of the method and reflects the required Hermite order, as described in the fol- 
lowing and thoroughly discussed in ref.— . The mesh is a Cartesian grid and the lattice 
velocities satisfy the sum rules, J2 p w p c pi = 0> 12 P w p c pi c pj = v< r^iji Yli P w p c pi c pj c pk — and 
^2pW p c P iC P jC p kCpi = v^(SikSji + SuSjk), in order to guarantee mass and momentum conserva- 
tion and isotropy. 

The third step of the procedure is to propagate the distribution via a discretization of 
the streaming operator to first order, as a forward Euler update, 

- 1 

f p (r + Atc p , t + At)= f p (r, t) + Atw p ]T — [^ l \r, t) + X «(r, *)] h^(c p ) (63) 



=o ° tL - 
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where At is the LB time-step. Consequently, the algorithm exploits the same Cartesian 
mesh to rearrange populations over spatial and velocity shifts. 

The population dynamics is able to reproduce the target macroscopic evolutions, such as 
the Navier-Stokes equation, to high accuracy. In the appendix a formal Chapman-Enskog 
multiscale analysis shows that for a second-order of accuracy in the transport coefficients, a 
second order expansion in Hermite is needed for both the equilibrium distribution and the 
collisional integral. As a result, the final form of the LB algorithm reads 

f p (r + c p , t + 1) = (1 - — )/ p (r, t) + — ft* (r, t) + Ai/Q(r, t) (64) 

T T 

where the relaxation time r is related to the kinetic component of the shear viscosity via 

rjW = n 4( r _ ^) (65) 

i.e. the physical value (compare with eq. (HDjl ) minus a viscosity of numerical origin. Some 
straightforward algebra shows that the final form of the discrete equilibrium reads 



fl q = w p n{r,t) 



c pi Uj(r, t) (c pi c pj - v^5ij)ui(T, t)uj(r, t) 

' 2 o 4 



(66) 



which is tantamount to a low-Mach (0[Ma 3 ]) expansion of the local Maxwellian, while the 
discretized form of the collisional term is 



/C p = -w p — 

m 



c pi C?\r,t) (c^-^^MjCfM) 

2 ' 4 



(67) 



A popular mesh model that reproduces the second order Hermite accuracy is provided by 
the so-called D3Q19 model^l, consisting of 19 discrete velocities in three dimensions, and 
with vt — 1/ V^- 

A further issue concerns the calculation of the spatial convolution present in the collisional 
term. From the above discussion it is clear that the central issue in LB related methods 
is the discretization of velocity space with the ensuing level of accuracy in the macroscopic 
transport equations. On the other hand, the method is completely flexible in terms of the 
mesh spacing Ax, that can be tuned at will in order to resolve the details of the microscopic 
interactions. Therefore, the error introduced in the spatial discretization does not represent 
a critical issue. As an example, in a previous paper 3 - we solved numerically eq. (1331) by 
means of the Lattice Boltzmann method 3 ^ using as a test case the flow of a hard-sphere 
fluid through a narrow channel with a pressure gradient along its axis, and measuring the 
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deviations from the macroscopic Poiseuille law. The spatial convolutions were evaluated with 
spatial quadratures over a number of off-mesh points obtained via linear interpolations, one 
of several alternatives to compute such integrals to desired level of accuracy. 

Regarding the stability of the proposed method, it is worth mentioning that even for 
simple non-interacting dynamics the LB algorithm is subjected to numerical instability. By 
using a von Neumann linear stability analysis, it has been shown that a stable LB scheme 
requires that the flow velocity be below a certain threshold that is a function of the relaxation 
time and the wave number—. The action of stiff intermolecular forces clearly narrows the 
stability range. Nevertheless, a generic upper bound for the variation of populations due to 
the forcing term is Sf/f ~ f P /w p ~ AtJC p <C 1. Therefore, for a generic quadrature scheme, 
the convolution force is C( l > / Ax 3 ~ nFg 2 , and the forcing term should be nFg 2 <C I/Arc 4 
since Ax/ At ~ 1, showing that the stability upper bound raises rapidly when reducing the 
mesh spacing. 

VI. CONCLUSIONS 

To summarize, we have presented a theoretical analysis and proposed a computational 
scheme which bridge hydrodynamics with microscopic structural theories of fluids. The 
present approach shows that the dynamic density functional and Boltzmann-like methods 
can be derived in a unique framework. The differences between the two methods are de- 
termined by the interaction of the fluid with the heat bath. The DDF method applies to 
colloids which are embedded in a solvent whose degrees of freedom are represented by a 
viscous heat bath, not moving with the particles, which eliminates all the hydrodynamic 
modes, but the diffusive mode. Boltzmann methods, instead, apply to molecular fluids and 
the corresponding heat bath is determined by degrees of freedom internal to the fluids and 
therefore moves with the fluid. 

A further improvement could concern the dependence of the hydrodynamic fields u(r, t) 
and T(r, t) on the local distribution /(r, v, t). One could expect that these fields should 
be obtained via a coarse-grained prescription starting from (/>m(t, v,t) because the hydro- 
dynamics cannot be extended beyond the size of few diameters. In fact, at scales shorter 
than a few molecular diameters hydrodynamic modes are strongly damped, with a decay 
time so short that the time autocorrelation function of a single molecule is isotropic. There- 
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fore, descriptions based on hydrodynamic field could become inappropriate at that scaled 
This program has been carried on for the hard sphere model in ref.— by employing a radial 
distribution function constructed according to the prescription of Fischer and Methfessel^. 
It is assumed that g 2 (r,r'\n) depends on a a coarse-grained density n(r,i) via a uniform 
smearing over a sphere of radius a/2. As a consequence the collisional terms (1461) and (14"T|) 
depend non locally not only from the density as in DDF but also from the temperature and 
velocity field. 

In the final section of the paper we showed how the derived kinetic equations can be 
transformed into a numerical scheme by following the theoretical derivation, i.e. by relying 
on a truncated Hermite representation of the collisional integrals and the unknown distri- 
bution, complemented by Gauss-Hermite quadratures to evaluate the distribution moments. 
Previous implementations of these ideas by us have proved that this numerical approach can 
be successfully applied to the study of inhomogeneous fluids. 
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APPENDIX A: CHAPMAN-ENSKOG ANALYSIS 

Let us analyze the multi-scale dynamics in powers of the Knudsen number e, with the 
distribution expanded as 



The dynamics is further decomposed based on the intrinsic hierarchy of timescales. For 
example, one can distinguish between the short-time convective motion, over a spatial scale 
x — > ex and temporal scale t± = et, and the slower evolution set by momentum diffusion, 
t 2 = e 2 t. Expressing the derivatives as 



/ = r ? + e/« + e 2 / {2) + ... 



(Al) 



d t = ed? + e 2 d? ) + 



(A2) 







eD 



(A3) 



the Lattice Boltzmann streaming step is expressed as 



f p ( r + Cp ,t+l)-f p (r,t) = (edl 1} +e 2 dP + ... + ec p d)(f; q + ^ ) + -) 
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having set At = 1. We now employ the explicit form of the kinetic equation (1331) . The 
collision operator has an intrinsic dependence on e, and the BGK term reads 

-<J?-fp) = + + (A4) 

T T 

where, as shown next, the relaxational time r is related to the kinetic component of the shear 
viscosity. The collisional term in (1331) is a non-linear functional of the density, multiplied by 
d v f eg , that depends on density and momentum density, so that the fastest evolution is on 
the e scale, 

/C = e/C« +e 2 K,M + ... (A5) 
and fC^ ' = 0. By equating the same orders in e, the evolutions on the e scale is 



{d[ x) + Cp d)f;« = + (A6) 

and on e 2 scale 



(i) 



T P 



dl 1] fW + df ] r v q + dcpf® + \c v c v d 2 r v q + ddl 1] c p f;« + ffiaPf? 

f (2) 

= -lZ- + lC® (A7) 
r F 

Substituting eq. (IA6l) into the fourth and last terms in the l.h.s. of eq. flA7j) . and by rearrang- 
ing terms, the e 2 dynamics can be rewritten as 

Up to the momentum diffusivity scale, therefore, the evolution of the density and current 
are obtained by contracting the above equations as J2 P (') an d J2 P c p(')y obtaining 

d^n + d^inui) = (A9) 
df ] n = (A10) 
d^inui) + dj (nv^Sij + n Ui Uj + P t f ,l) ) = (All) 

d? (n Ui ) + \ dj ( (l - ^) if + 4 C ' 2) ) = (A12) 

where P^? and P^' 2 ^ are the collisional contributions to the pressure tensor at first and 
second Knudsen order, respectively In addition, some lengthy algebra shows that p^ ,neqS> 
is identified as the non-equilibrium part of the kinetic pressure tensor 

Pf neq) = (l - ^) E - ^ = V {K) [diUj + d jUi ] (A13) 
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By multiplying eqs. (|A9]) and flAlOH by e and eqs. (lAlip and (IA12D by e 2 and summing the 
equations, we reconstruct the sought evolution for density and momentum density accurate 
to e 2 level 



d t n + dAnui 



d t {n Ui ) + d j (nv 2 T 6 tJ + Pl K) + P l 



(C) 

ij 









At higher Knudsen order, the following recurrence relation holds 



m=0 



(A14) 



where = / e9 . Similarly, for the Hermite coefficients, 



(n,fc+l) 



E 

m=0 



(A15) 



where the first superscript of the coefficients refers to the Hermite truncation level and the 
second to the Knudsen level. If we neglect the collisional term, eq. (1A15|) shows that the 
dynamics of the coefficient 0( n ' fc+1 ) depends on cf)^ n ~ 1 ' k \ 0( n+1 > fc ) an d (p^ n '°\ 0( n,fc ) in a 
hierarchical (pyramidal) way. In order to guarantee a given level of accuracy on the Knud- 
sen level, without introducing any approximation from the Hermite truncation, a complete 
representation is required on the base of the pyramid corresponding to the equilibrium dis- 
tribution. In particular, for Knudsen order k, the equilibrium needs to include at least 
n + k + 1 Hermite coefficients 42 . At first order in Knudsen, the equilibrium distribution 
should be accurate up to third Hermite order, complemented by fifth order quadratures to 
evaluate the moments. In order to handle thermal dynamics (in the presence of a non-zero 
heat flux), a fourth order accuracy in the equilibrium is required, together with seventh 
order quadratures. 

When the collisional term is included, the simple hierarchical picture seems to be spoiled. 
However, by using the explicit Hermite representation of the collisional term, the third order 
Hermite coefficient ~ u 3 , a contribution that can be neglected in standard condensed 
matter conditions. If one retains only up to second order in fp eq \ the error can be estimated, 
as discussed irA to be ~ —u 2 d x u in the pressure tensor and ~ u 2 in the shear viscosity. 
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